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LATTICE BOLTZMANN EQUATION ON A 2D RECTANGULAR GRID 

M ’HA MED BOUZIDP, DOMINIQUE D’HUMIERESt, PIERRE LALLEMAND*, AND LI-SHI LUO$ 

Abstract. We construct a multi-relaxation lattice Boltzmann model on a two-dimensional rectangular 
grid. The model is partly inspired by a previous work of Koelman to construct a lattice BGK model on a 
two-dimensional rectangular grid. The linearized dispersion equation is analyzed to obtain the constraints 
on the isotropy of the transport coefficients and Galilean invariance for various wave propagations in the 
model. The linear stability of the model is also studied. The model is numerically tested for three cases: (a) 
a vortex moving with a constant velocity on a mesh with periodic boundary conditions; (b) Poiseuille flow 
with an arbitrary inclined angle with respect to the lattice orientation; and (c) a cylinder asymmetrically 
placed in a channel. The numerical results of these tests are compared with either analytic solutions or the 
results obtained by other methods. Satisfactory results are obtained for the numerical simulations. 

Key words, generalized lattice Boltzmann equation, rectangular meshes, stability analysis, dispersion 
equation, Taylor vortex, Poiseuille flow, flow past a cylinder in channel 

Subject classification. Fluid Mechanics 

1. Introduction. Historically originating from the lattice gas automata (LG A) introduced by Frisch, 
Hasslacher, and Pomeau [3], the lattice Boltzmann equation (LBE) has recently become an alternative 
method for computational fluid dynamics. The essential ingredients in any lattice Boltzmann models which 
are required to be completely specified are: (i) a discrete lattice space on which fluid particles reside; (ii) a 
set of discrete velocities (often going from one node to its nearest neighbors) to represent particle advection; 
and (iii) a set of rules for the redistribution of particles residing on a node to mimic collision processes in 
a real fluid. Fluid-boundary interactions are usually approximated by simple reflections of the particles by 
solid interfaces. 

In a hydrodynamic simulation by using the lattice Boltzmann equation, one solves the evolution equations 
of the distribution functions of fictitious fluid particles colliding and moving synchronously on a highly 
symmetric lattice space. The highly symmetric lattice space is a result of the discretization of particle 
velocity space and the condition for synchronous motions. That is, the discretizations of time and particle 
phase space are coherently coupled together. This makes the evolution of the lattice Boltzmann equation 
very simple — it consists in only two steps: collision and advection. One immediate limitation of the 
LBE method is due to its use of highly symmetric regular lattice mesh, which are usually triangular or 
square lattices in two dimensions and cubic in three dimensions. Obviously this is a serious obstacle to its 
applications in many areas of computational fluid dynamics. To deal with complex computational domains, 
various proposals have been made to use grids that are better suited to fit boundaries or to adapt meshes 
according to the physics of the system. 
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It has been shown recently that the lattice Boltzmann equation is indeed a special finite difference 
form of the continuous Boltzmann equation with some drastic approximations tailored for hydrodynamic 
simulations [4, 5, 7]. This makes the lattice Boltzmann method more amenable to incorporate body- fitted 
meshes [8, 9] or grid refinement techniques [10]. In most cases, the regular lattice mesh is abandoned by 
decoupling the spatial-temporal discretizations and the discrete velocity set. so that interpolations can be 
used in addition to the advection on a non-regular or non-uniform mesh. However, interpolations introduce 
additional numerical viscosities and other artifacts into the lattice Boltzmann method [11]. Therefore, it is 
highly desirable to construct lattice Boltzmann models with arbitrary mesh and free of interpolations [2, 12]. 

In this paper we shall consider a two dimensional model on a rectangular grid with an aspect ratio 
of a — Sy/&x , where 0 < a < 1. The model is inspired in part by a previous work of Koelman [2] who 
proposed a general scheme to construct lattice BGK models with given discrete velocity sets based on a low 
Mach number expansion of the Maxwellian equilibrium distribution function. Conservation and symmetry 
constraints are imposed to fix the parameters in the equilibrium distribution function. Koelman ’s model is 
essentially a variation of the lattice BGK model [13, 14]. As we shall show, the transport coefficients of this 
model are generally anisotropic when a ^ 1 [15]. 

We use the generalized lattice Boltzmann equation with multiple relaxation times due to d'Humieres [1], 
instead of the standard lattice BGK model [13, 14]. The generalized LBE model has the freedom of multiple 
relaxations which can be independent or coupled together. This allows one to optimize the overall properties 
of the model through suitable compensation of inadequate behaviors. We shall study the time evolution 
of plane waves by analyzing the linearized dispersion equation of the model [11]. This analysis allows us 
to obtain the conditions under which the model can be used to simulate the Navier-Stokes equation, i.e., 
the model is Galilean invariant and isotropic up to a certain order in wave-number k. We show that severe 
stability constraints are due to Galilean invariance and isotropy of transport coefficients. This demonstrates 
the difficulty in the endeavor of constructing a lattice Boltzmann model with arbitrary grid. Simulations of 
non-trivial cases are presented to demonstrate the qualities and defects of the model. 

We organize the paper as follows. Section 2 describes the proposed model on a rectangular grid. Section 3 
shows a detailed analysis of the dispersion equation. The wave-number dependence of Galilean coefficient and 
attenuation coefficients are computed explicitly to obtain the conditions under which the model is Galilean 
invariant and isotropic. Section 4 provides examples of numerical simulations using the proposed model: (a) 
a vortex moving with a uniform velocity in a periodic system; (b) Poiseuille flow with the boundaries along 
arbitrary direction with respect to the underlying lattice; and (c) flow past a cylinder asymmetrically placed 
in a channel. Section 5 concludes the paper. 

2. Definition of the model. We consider a two-dimensional LBE model with nine discrete velocities 
(the D2Q9 model) on a rectangular grid of dimensions 1 and a. (In what follows all quantities are given in 
non-dimensional units, normalized by the lattice unit <5*.) In the advection step of the lattice Boltzmann 
equation, particles move from one node of the grid to one of its neighbors as illustrated in Fig. 1. The 
discrete velocities are given by 


e Q 


r (0, 0) , q = 0 , 

< (cos[(q - l)7r/2], a sin[(a - 1 )tt/2]), a = 1 -4, 

_ (cos[(2a — 9)7r/4], a sin[(2a — 9)7 t/ 4])\/2, a = 5-8, 
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where the duration of the time step St is assumed to be unity. At any time t n , the LBE fluid is then 
characterized by the populations of the nine velocities at each node of the computational domain 

I f(rj, t n )) = (fo(rj, «„), t n ), ••• , t n )) T , (2.2) 

where T is the transpose operator. Here upon the Dirac notation of bra, (*|, and ket, |), vectors is used to 
denote row column and row vectors, respectively. The time evolution of the state of the fluid follows the 
general equation 

\f(rj + e a . t n + 1)) = |/(r i , <„)) + \n(f(r jy t n ))} , (2.3) 


where collisions are symbolically represented by the operator 0. 

h 2 



e 7 e 4 e 8 


Fig. 1. Discrete velocities of the nine-velocity model on a rectangular grid . 


The aspect ratio of the rectangle is Syj&x = a. 


We shall use the generalized lattice Boltzmann equation introduced by d’Humieres [1], in which the 
collision process is executed in moment space M. The mapping between moment space and discrete velocity 
space V spanned by {e a } is one-to-one and defined by the linear transformation M which maps a vector |/) 
in V to a vector |/) in M, he., 

|/) = M|/), and |/) = M- 1 |/). (2.4) 


To reflect the underlying symmetries appearing in both the Chapman-Enskog expansion and the dispersion 
equation, M is constructed as the following 
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= (|mi), | m 2 ), | m 3 ), |m 4 ), |m 5 ), |m 6 ), |m 7 ), |m 8 ), |m 9 )) T , 


(2.5) 


where = a 2 + 1, <^2 = 1 — 2a 2 , <£3 = a 2 — 2, <£4 = a 2 — 1, ip $ = a 2 + 2, and — —(1 + 2a 2 ). 

The components of the row vector (mp\ in matrix M are polynomials of the x and y components of the 
velocities {e Q }, e a , x and e a , y . The vectors (mp |, 0 = 1, 2, * ■ • , 9, are orthogonalized by the Gram-Schmidt 
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procedure in a carefully considered order. The first three orthogonal vectors correspond to the mass, x- and 
^-momentum modes: (mi| = (||c Q ||°|, (m 4 ] = (e a . x |, and (m 6 | — (e Q>y |. The above expressions prescribe the 
components of (mi |, {m±\, and (me |. These three vectors span the hydrodynamic subspace of the eigen-space 
of the collision operator for a two-dimensional athermal LBE model. The remaining six vectors span the 
kinetic subspace. The vector {m 2 | = (3||e a || 2 — 2(1 -f a 2 )||e a ||°| is constructed by orthogonalizing the energy 
mode (||e a || 2 |. Similarly, vectors (m 5 | and {m^\ are respectively built upon (e Q?x ||e a || 2 | and {e Q , y ||e Q || 2 |; 
(mg | is constructed upon (e 2 x — e 2 y | and (mg | = (e a?x e Qiy |; and finally (m 3 | is obtained from (||e a || 4 |. By 
means of their construction, the row vectors in M are mutually orthogonal, but they are not normalized, 
their norms being chosen to simplify algebraic manipulations. When a .= 1 , M reduces to that for the D2Q9 
model on a square grid with a different normalization of | p xx ) [11]. Therefore the proposed model can be 
considered as a generalization of the model on a square lattice. It should be noted that when a ^ 1 , there 
are three nonzero (kinetic) energy levels in the model which introduce additional degrees of freedom into 
the model and extra care must be taken in the construction of (m 2 |, (m 8 |, and (m 3 |, i.e., they must be 
orthogonalized with the Gram-Schmidt procedure in the particular order of (ra 2 |, (mg|, and (77i 3 |. 

It is interesting to note that the moments (m$ \f) have a physical interpretation. The matrix M so 
constructed in the above naturally leads to the moment vector in moment space M as the following: 

I/) = (.Pt jx 1 Qxy jy 5 Qyy Pxx) Pxy ) •> (2-6) 


where p is the density, e is related to the kinetic energy, e is related to the kinetic energy squared for a — 1 
(but has no obvious physical meaning when a / 1 ), j x and j y are x and y components of the momentum 
density, q x and q y are proportional to the x and y components of the energy flux, and p xx and p xy are 
proportional to the diagonal and off-diagonal components of the viscous stress tensor. 

For the collision process, we propose to use the following equilibrium distribution functions of the 
(nonconserved) moments, which depend only on the conserved moments, i.e., p, and j = ( j x , j y ): 


5 < e< i> = 2(3Cj — l — a?)p + -{jl + jl) , 


1 


£ <eq) = ;| q 3P, 

?i eq) = Tux, 


% 


(eq) _ 1 


C 2 Jy 


ifi?’ = [ 3(« 2 + D4 - 2« 2 ] p + l (fjl - ij 2 ) , 


P*v ] ~ pixiy ■ 


(2.7a) 

(2.7b) 

(2.7c) 

(2.7d) 

(2.7e) 

(2-7f) 


where the coupling coefficient between p ^ and p (which vanishes in the standard D2Q9 LBE model) is 
introduced to obtain the isotropy of the sound speed. The values of the coupling constants (a 3 , c\, and c 2 ) 
in the above equilibria are obtained by optimizing isotropy and and stability of the model [11]. It should be 
noted that the energy is not considered as a conserved quantity here because the model is athermal. (The 
model does not possess sufficient degrees of freedom to accommodate the dynamics of locally isotropic heat 
transport). 

In what follows the idea of the “incompressible” LBE [ 6 ] is applied to the above equilibria so that p 
is replaced by a constant po in the denominators of equations (2.7a), (2.7e), and (2.7f). This choice allows 
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for better comparison with other incompressible simulations and simpler algebra while retaining correct 
acoustics. 

The collision process is modeled by the following relaxation equations 

in = i/Vs[i/Vi/ (eq, )]i ( 2 - 8 ) 

where |/*) denotes the post-collision state, and S is the diagonal relaxation matrix 

S = diag(0, $ 2 , 0, s 5 , 0, s 7 , s s , s 9 ) . (2.9) 

The model reduces to the usual lattice BGK model if all the relaxation parameters are set to be a single 

relaxation time r (and a - 1), i.e., s a = 1/r. It should be stressed that the relaxation parameters are not 
independent, as shown in the next section. The constraints of isotropy lead to the coupling between these 
relaxation parameters [11]. Obviously, the usual lattice BGK model does not possess the freedom for such 
couplings, therefore it would not work on a rectangular grid. 

3. Analysis of linearized dispersion equation. The analysis presented in what follows is similar 

to that presented in Ref. [11] where the goal of the work was to determine the stability conditions for the 

coupling coefficients a 2 and e* 3 , and the constraints on the relaxation parameters s a . 

We consider a system of size N x x N y with periodic boundary conditions and look for small amplitude 
solutions in the presence of a uniform flow [for given values of p and V = (V x , V y ) = J / p\. For a wave vector 
k in the reciprocal space of the computational domain, we search for solutions 

f a (r, t) ex exp(-zfc • r + zi) . (3.1) 


To first order in fc, we have the following linearized dispersion equation: 

det(K (1 ^ + M -1 CM - z\) — 0, 


(3.2) 


where I is the identity operator, is the linearized advection operator which is a diagonal matrix 


K (1) = diag(0, ik * e \ , . . . , ik • eg) , 


and C is the linearized collision operator 
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Q 4 = 2[3 Cg - (1 + a 2 )] , 

a s = ^^[2a 2 -3cV + l)] . 
a 1 


(3.3) 


(3.4) 


(3.5a) 

(3.5b) 
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The linearized dispersion equation (3.2) can be solved by perturbation technique in power series of k [11]. 
To ensure isotropy and Galilean invariance in the limit of k — ► 0, we need to solve the linearized dispersion 
equation up to k 2 . 

In the first order of k , three solutions are obtained: one corresponds to transverse excitations which are 
convected with the uniform speed of the fluid k • V /k y whereas the other two are acoustic waves with phase 
velocity ±c 8 , where the speed of sound c s can be chosen within limits that is deferred to later discussion. 
The sound waves also have the correct dependence on the applied uniform velocity V of the fluid up to the 
first order in V, i.e., c 8 -4 c s ± Vcos</>, where <p is the angle between k and V . The nonlinear terms in the 
equilibria of Eqs. (2.7a) -- (2.7f) provide the correct Galilean coefficients for both transverse and longitudinal 
waves. 

In the second order (in k) of the solutions of the dispersion equation, the constraints on the isotropy of 
the transport coefficients for the hydrodynamic modes lead to 


c 2 + 4(1 - a 2 ) 

Cl = 1 

a 1 

and the following relationships between the relaxation parameters 

1 _ 2(4+c 2 )[(12 C 2- C2 )(l+a 2 )-2(5a 2 + 2)] 1 

s 2 (l + a 2 )(l + c 2 — 3a 2 )(c 2 -j- 10— 12c 2 ) + 6[a 4 (c 2 — 2) — 3(a 2 — 1)] s 9 5 
1 2(4+c 2 )[(12c 2 -c 2 )(l-fa 2 )-2(3a 4 -h5a 2 -h5)] 1 

s 8 “ (l+a 2 )(l + c 2 — 3a 2 ) (c 2 + 10— 12c 2 ) +6[a 4 (c 2 — 2)— 3(a 2 - 1)] I9 ’ 


(3.6) 


(3.7a) 

(3.7b) 


where 1 fs a = (1 /s Q - 1/2). The coupling between s 2 and s 9 is required only when a ^ 1. The kinematic 
shear viscosity u and the kinematic bulk viscosity £ are 


4 + c 2 (1 1\ 

6 Us V ’ 


<=^(r + 3a ! + c 2 -12c;) (i-i) ■ 


(3.8a) 

(3.8b) 


For a given a, the speed of sound and c 2 must be chosen such that the shear and bulk viscosities are positive 
and the Eqs. (3.7a) and (3.7b) lead to positive values for s 2 and s 8 . 

The values of c 8 and c 2 which optimize the isotropy and stability of the model, depending on the grid 
aspect ratio a, are determined by the linear analysis of the model [11]. In the case of square grid, i.e., 
a = 1, we have found c 2 = 1/3 and c 2 = —2. This result coincides with the one given in Ref. [11] and the 
relationship between and s 9 given by Eq. (3.7b), and the shear and bulk viscosities given by Eqs. (3.8a) 
and (3.8b), all reduces to the previous results for a square lattice where c\ = — 2 [see Eqs. (40), (41), (42), 
and (43) in Ref. [11] for c s , sg(s 8 ), v , and £, respectively]. However, the coupling between s 2 and s 9 is unique 
to the model on a rectangular grid. This coupling is due to the dependence of p ^ on p, which in turn leads 
the term a 5 s 8 in the linearized collision operator C in Eq. (3.4). Finally note that q 3 has little influence and 
is set to be equal to —2. 

The linearized dispersion equation can be solved numerically for any value of k to determine the linear 
stability of the system by computing the rate of growth of spatially periodic excitations superimposed to 
a uniform flow of constant velocity V, as previously shown in the case of a square grid [11]. Through this 
analysis it is found that the present model is much less stable than the square one, i.e., the stable region in 
parameter space of V and s a is much smaller than that for the model with a square lattice. For instance, 
when a = 1/2, a stability condition is that V < 0.05, whereas for the model with a square lattice (a = 1), the 
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same stability condition is that V < 0.20. One reason for this is that in general the sound speed c s decreases 
with the aspect ratio a; for instance, when a = 1/2 the optimal speed of sound is about 0.377. which is 
different from the usual c s — l/\/3 % 0.577 on a square lattice. Therefore, the local velocity magnitude 
must be decreased accordingly to keep the local Mach in check so that the low Mach number approximation 
remains valid. This means that the present model will have limited ability to simulate flows even at moderate 
Reynolds numbers. In addition when using a combination of rectangular and square grids (the simplest case 
of grid refinement in one direction) in the situation where acoustic propagation is important, it will not be 
possible to choose an optimal value of the sound speed for the two different grids. 

We would like to note that, although there is no simple interpretation of the instability of the LBE 
models due to the presence of a uniform velocity V, information on the instability can be obtained by 
analyzing the velocity dependence of the attenuation of sound waves using the linearized dispersion equation 

in]. 

Let us consider the case where the uniform velocity is parallel to the wave vector k with a polar angle 
6 (between k and ar-axis). For small values of k and the particular choice of co 

c 2 = (a 2 - 3), (3.9) 


we have the following results. The transverse mode has phase velocity v± = V and its attenuation is given 
by 


7-L 




9(1 - a 2 ) 2 sin 2 28 
2[1 - 13a 2 + a 4 +6(1 + a 2 )cf] 



(3.10) 


For the longitudinal modes, we obtain as phase velocity uy = ± + V 2 and attenuation coefficient 7|| = 
{lb + 7±)/2, with (to the first order in V r ) 




3c: 


V 


3 c s [l+8a 2 + a 4 -12(l+a 2 ) c 2j 

-3 (7 + 12a 2 + 3a 4 ) c 2 +12(1 - a 2 )c 2 cos 2 0 [2 + (1 - a 2 )(2 - 3cos 2 (9)]} ) . 


{(1 + a 2 )(7a 2 + 36c 4 ) 


(3.11) 


Contrary to the case of square grid, it is not possible in general with a given value of a ^ 1 to find a value of 
c s for which the linear dependence of the attenuation of acoustic waves on V can be eliminated (for a = 1, 
this can be accomplished by setting c 2 = 1/3). This is a possible cause of instability in the model. 

4. Simulations. We use the two-dimensional multi-relaxation LBE model on either a square grid or 
a rectangular grid for the following simulations. The central routine (collision and advection) is quite close 
to that for the standard square LBE and leads to similar performances (using a workstation with a 500 
MHz EV6 processor, the overall computation time per node and per time step is in the range 0.2 to 0.4 
microsecond depending whether the cache is large enough or not). 


4.1. a vortex traveling with a constant velocity. To test the ability of the present LBE scheme to 
simulate a viscous flow, we consider the particular case of a simple vortex superimposed to a uniform flow 
of velocity V. We take as initial condition for the flow 


u 0 (r, t = 0) = V + (yo - y, x -i 0 )w 0 exp[-(r - r 0 ) 2 /ft 2 ]. 


(4.1) 


where ro = (:ro, Vo) is the initial position of the vortex center, and cuo and R characterize respectively the 
amplitude and the extent of the vortex. The evolution of the corresponding macroscopic flow is fairly simple: 
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the center of the vortex travels with the velocity V and the maximum value of the vorticity (at the center 
ro + Vt) decays in time as 

= (F + 4°^) 2 = (1+4 i*) 2 ’ (4 ' 2) 

where t* = ut/R 2 is the dimensionless time. 



Fig. 2. LBE simulation of a moving vortex. Decay of the vorticity maximum. The grid aspect ration a ~ 1/2. Symbol + 
and x are simulation results for V — 0 and V = 0.05, respectively. The solid lines are fitting of the data according to Eq. (4-2) 
with the viscosity of value v — 0. 9876^0 and v — 0. 8966^0, respectively, where v o is given by Eq . (3.8a). 

The system size is N x x N y = 109 x 109, with a grid aspect ratio a = 1/2. The size of the vortex is 
R — 6. Values of other parameters are: Q 2 = —3.5, a 3 = 2.0, C 2 = —2.9, and 58 = 1.8, Le., v = 0.01018 
according to Eq. (3.8a). The results obtained by the LBE simulations with various conditions agree very well 
with the analytic solution of the flow for V = 0. However when V increases there are departures from the 
simple result of Eq. (4.2) due to the dependence of the transport coefficients and ^-factor on V, as discussed 
for the square grid in Ref. [11]. An example of such behavior is demonstrated in Fig. 2. Figure 2 shows 
two LBE simulation results of u; ma x as a function of dimensionless time t* = ut/R 2 , with V = (0, 0) and 
V = (0.05, 0). Equation (4.2) is used to fit the data to obtain the viscosity. The results are v = 0.9876t/o 
and v — 0.8966^o for V x = 0 and V x = 0.05, respectively, where i/ 0 is given by Eq. (3.8a). There are two 
factors that attribute to the correction in the viscosity: the wave-number fc-dependence and V-dependence 
of the transport coefficients [11]. The same simulations are performed on a square grid and the results are: 
v = 0.9866j/o and v — 0. 8745^0 for V x = 0 and V x = 0.05, respectively. It should be noted that in the LBE 
simulations, initial conditions include not only the conserved quantities such as the density and velocity 
fields, but also all the nonconserved quantities such as fluxes and the stress, which can be obtained from the 
initial velocity field through a Chapman-Enskog analysis of the model. 

4.2. Poiseuille flow with arbitrarily inclined walls. The second test is the two-dimensional Poiseuille 
flow with arbitrarily inclined walls. This situation allows us to test the no-slip boundary conditions in the 
LBE model. We consider a system of size N x x N y with periodic boundary conditions. The boundaries of 
the channel are placed with an arbitrary inclined angle 0 with respect to z-axis, as illustrated in Fig. 3. 
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The no-slip boundary conditions used here for the channel walls are the interpolated bounce-back boundary 
conditions proposed in Ref. [16]. The interpolated bounce-back boundary conditions combine interpolation 
and bounce-back schemes to deal with boundaries which are off the lattice points. 



Fig. 3. 2D Poiseuille flow with arbitrary inclined walls. The system size is assumed to be N x x Ay. The discs are grtd 
points. The solid lines are the advection lines of the discrete velocities . The dashed lines are the boundaries of the channel. 
The width of the channel is N y . The no-slip boundary conditions are enforced at the intersections the dashed lines and the 
thin solid lines. 

We first studied the time evolution of the flow starting at rest, and compared the results obtained by 
using the rectangular and square grids. The time evolution of velocity fields of the two systems agree very well 
with each other. We also studied the momentum transfer at the boundary. We found an excellent agreement 
between its measurements for the square and the rectangular grids, and its expected value: pvLd j_ Vjj , where 
L is the length of the boundary, and d±V\\ is the normal derivative of the shear velocity with respect to the 
wall, computed at the wall. 

Note that when we compute the momentum transfer for the rectangular grid, the components of the 
usual momentum transfer have to be multiplied by a factor a to account for the surface of the elementary 
cell (assuming that all results are in non-dimensional units defined on the square grid). In order to better 
understand the origin of this factor, one has to remember that p and j are the mass and momentum densities 
(mass and momentum per unit surface), while the momentum transfer has to be computed from momentum: 
(momentum density) x (cell surface). Usually a unique regular grid is used and the cell volume can be taken 
as unit volume. Here however the surface of the cells is equal to the chosen aspect ratio a once the square 
cell has been taken as unit surface. Indeed this remark applies to the next section when computing the drag 
and lift coefficients. 

4.3. flow past a cylinder asymmetrically placed in a channel. The third test we did was a two- 
dimensional flow past a cylinder asymmetrically placed in a channel. This flow T has been used as a standard 
benchmark test in CFD [17]. The flow configuration is as follows: a cylinder of diameter d is placed in a 
channel of width 4. Id and length 22 d, the center of the cylinder is asymmetrically (with respect to the center 
line of the channel) located at horizontally 2d from the entrance, and vertically 2d from the lower wall of the 
channel, as shown in Fig. 4. The boundary condition at the entrance is a Poiseuille profile with average speed 
U. The boundary condition at the exit is free exit with a total flux equal to the input flux. The bounce-back 
boundary conditions are used for the channel walls, and the interpolated bounce-back boundary conditions 
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Table 4.1 

2D flow past a cylinder asymmetrically placed in a channel at Re — 100. 


a 

C 8 

C2 

Ax X Ny 

St 

fmax 

r* max 

l l 

/'“nil in 
L L 

AP 

1.00 

1/V3 

-2 

709 x 132 

0.3021 

3.153 

0.926 

-1.018 

2.50 

0.85 

0.6141 

-0.80 

709 x 155 

0.3018 

3.186 

0.984 

-1.051 

2.51 

0.80 

0.5829 

-0.90 

709 x 165 

0.3020 

3.174 

0.950 

-1.062 

2.51 

0.75 

0.5412 

-1.10 

709 x 176 

0.3018 

3.173 

0.965 

-1.053 

2.51 

0.70 

0.5113 

-1.50 

709 x 188 

0.3007 

3.195 

1.013 

-1.071 

2.51 

0.65 

0.4761 

-1.70 

709 x 203 

0.3009 

3.184 

0.999 

-1.062 

2.47 

0.60 

0.4417 

-2.00 

709 x 220 

0.3009 

3.176 

1.002 

-1.053 

2.45 

0.55 

0.4086 

-2.25 

709 x 240 

0.3015 

3.189 

1.005 

-1.052 

2.42 

0.50 

0.3770 

-2.55 

709 x 264 

0.3007 

3.199 

1.019 

-1.084 

2.45 

0.45 

0.2977 

-2.90 

709 x 293 

0.2992 

3.204 

1.053 

-1.107 

2.50 

CFD lower bound in Ref. [17] 

0.2950 

3.22 

0.99 

— 

2.46 

CFD 

upper bound in ] 

Ref. [17] 

0.3050 

3.24 

1.01 

— 

2.50 


with a second order interpolation [16] are used for the boundary of the cylinder. The Reynolds number for 
the flow is 



v 


We use the LBE model to simulate the flow at Re = 100 for which there is periodic vortex shedding behind 
the cylinder. 




16 d 


d ) 

1 . 5d 


22 d 


(u. v) = <0, 0) 


flow direction ► 


(u, V ) = (O, 0) 


Fig. 4. Configuration of a 2D flow past a cylinder asymmetrically placed in a channel. 

The flow was computed on rectangular grids with several different values of the grid aspect ratio a, and 
compared to the results with a square grid. The measured quantities are Strouhal number St, maximum drag 
Cg ax , maximum lift coefficient minimum lift coefficient C™ m , and the pressure difference A P. The 

results are summarized in Table 4.1. Table 4.1 also shows the lower and upper bounds of St, Cj5 ax , , and 
A P, obtained by a number of conventional CFD methods presented in Ref. [17]. Overall the LBE simulation 
results with square or rectangular grids agree well with each other, and with the CFD results in Ref. [17]. 
Figures 5 show 7 the contours of the stream function ifcix, y) and the vorticity a>(;c, y) of the simulations on 
a square grid of size N x x N y = 1401 x 308 and on a rectangular grid of size N x x N y = 1401 x 616. The 
relative L 2 -norm difference of the two velocity fields is about 2.2 x 10 -4 . Note that the aspect ratio for this 
particular calculation is slightly different from that showm in figure 4, but this has negligible effect for the 
present purpose of comparing results on the square and the rectangular grids. 

The relative fluctuation of Strouhal number St is well under 1% and the values of St are well within the 
bounds in Ref. [17]. The fluctuation of C™ ax is also under 1% but the values of Cg ax are all slightly lower 
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than the results in Ref. [17]. The fluctuation of A P is about 1% and the values of A P agree well with the 
results in Ref. [17]. The values of lift coefficient obtained by the LBE simulations have a variation about 
±6%, which is much greater than the variations in other measured quantities. 



0 200 400 600 800 1000 1200 


Fig. 5. 2D flow past a cylinder asymmetrically placed in a channel at Re = 100. Top and bottom figure show contours of 
the stream function i j){x, y ) and the vorticity uj(x^ y) of the flow, respectively . The dashed lines are the simulation results on 
a square grid of size N x x N y = 1401 x 308, and the solid lines are that on a rectangular grid of size N x x N y = 1401 x 616. 

A possible origin of the discrepancy in the lift coefficients is the following. The LBE method is intrin- 
sically a compressible scheme and acoustic w r aves may be generated by, e.g., initial conditions that do not 
include a proper pressure field or the flow itself that generates an oscillating pressure field as is the case 
considered here. For a given value of the sound speed and a given choice of the boundary conditions at the 
entrance and exit of the channel the frequency of some of the longitudinal acoustic modes can be close to 
multiples of the Strouhal frequency in the flow. This causes resonances between some of the acoustic waves 
and the periodic shedding of vortices by the cylinder. The coupling between acoustic w ? aves and vortex shed- 
ding indeed affects the hydrodynamic fields, and in turn, various measured quantities. Among the measured 
quantities, the lift coefficients are most sensitive to this effect. The mean drag coefficient is also affected but 
to a much smaller extent. This problem is of broad interest. However it will be easier to study it w r ith the 
model of square grid for which the speed of sound and the bulk viscosity can be chosen in a broader range 
than for the model of rectangular grid. A detailed study is beyond the scope of the present work and wall 
be addressed elsewhere. 

5. Conclusion and discussion. In this paper we have successfully proposed a tw'o-dimensional nine- 
velocity generalized lattice Boltzmann model with multiple relaxations on a rectangular grid with arbitrary 
aspect ratio a = S y /S x . We have numerically validated the model by using the model to simulate several 
benchmark problems, and have obtained satisfactory results. In contrast to the previous two-dimensional, 
nine- velocity, multi-relaxation model on a square grid [11], the model on a rectangular grid is more prone 
to instability, and the admissible maximum value of local velocity magnitude is much less than that in 
the model on a square grid. It should also be stressed that, although this work is in part motivated by a 
previous work [2], it is realized that the nine- velocity lattice BGK equation cannot possibly work properly 
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on a rectangular grid. Specifically, the lattice BGK equation does not have sufficient degrees of freedom 
to satisfy the constraints imposed by isotropy and Galilean invariance. With nine discrete velocities in 
two-dimensions, it is necessary to use the multi-relaxations to construct an LBE model on a rectangular 
grid. 

This work is our first attempt to construct a lattice Boltzmann model on an arbitrary unstructured grid. 
As discussed in Ref. [12], one difficulty encountered in the LBE model on an unstructured grid is due to 
the fact that Ve a f e a V/ because the discrete velocity set {e a } has spatial dependence. In this work, we 
found that there are additional issues in the LBE model on an unstructured grid needed to be addressed. 

First, we found that the local grid structure severely affects the local sound speed. If the sound speed 
varies spatially depending on local grid structure, then the model is unphysical. Correct acoustic propagation 
is an essential part of the lattice Boltzmann method. Secondly, the constraints of isotropy and Galilean 
invariance are difficult to satisfy by using the lattice BGK model, as proposed in Ref. [12], unless the discrete 
velocity set includes a large number of velocities. Thirdly, the numerical stability is severely affected by the 
local grid structure even for uniform structured grid, as we have demonstrated in this work. Stability is of 
key importance to an effective lattice Boltzmann algorithm. However, we have not yet developed a method 
to systematically improve the stability of the lattice Boltzmann method. We believe that the aforementioned 
issues must be resolved before we can construct a lattice Boltzmann model on an arbitrary unstructured 
grid. 
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